Climate heterogeneity shapes phylogeographic pattern of Hippophae gyantsensis (Elaeagnaceae) in the east Himalaya‐Hengduan Mountains

Abstract The interaction of recent orographic uplift and climate heterogeneity acted as a key role in the East Himalaya‐Hengduan Mountains (EHHM) has been reported in many studies. However, how exactly the interaction promotes clade diversification remains poorly understood. In this study, we both used the chloroplast trnT‐trnF region and 11 nuclear microsatellite loci to investigate the phylogeographic structure and population dynamics of Hippophae gyantsensis and estimate what role geological barriers or ecological factors play in the spatial genetic structure. The results showed that this species had a strong east–west phylogeographic structure, with several mixed populations identified from microsatellite data in central location. The intraspecies divergence time was estimated to be about 3.59 Ma, corresponding well with the recent uplift of the Tibetan Plateau. Between the two lineages, there was significant climatic differentiation without geographic barriers. High consistency between lineage divergence, climatic heterogeneity, and Qingzang Movement demonstrated that climatic heterogeneity but not geographic isolation drives the divergence of H. gyantsensis, and the recent regional uplift of the QTP, as the Himalayas, creates heterogeneous climates by affecting the flow of the Indian monsoon. The east group of H. gyantsensis experienced population expansion c. 0.12 Ma, closely associated with the last interglacial interval. Subsequently, a genetic admixture event between east and west groups happened at 26.90 ka, a period corresponding to the warm inter‐glaciation again. These findings highlight the importance of the Quaternary climatic fluctuations in the recent evolutionary history of H. gyantsensis. Our study will improve the understanding of the history and mechanisms of biodiversity accumulation in the EHHM region.


| INTRODUC TI ON
The East Himalaya-Hengduan Mountains (EHHM) is one of the most well-known biodiversity hotspots in the world (Myers et al., 2000).
This region is not only the diversification or/and origin center of many taxa (Hu, 1994;Jia et al., 2012;The Comprehensive Scientific Expedition to the Qinghai-Xizang Plateau, 1986;Wu, 1988) but also has been considered a major glacial refuge during the glacial period for many species (Frenzel et al., 2003;Wu, 1988;Yang et al., 2008).
Over the last two decades, this region has become recognized as an important issue in biodiversity and evolutionary research (Favre et al., 2015;Qiu et al., 2011;Renner et al., 2016;Spicer et al., 2021;Wang et al., 2011;Wen et al., 2014). The EHHM region has extremely complex topography and variable climate, and climatic gradients are steep from west to east (The Comprehensive Scientific Expedition to the Qinghai-Xizang Plateau, 1997). Also, here orogenesis has been strong since the early Cenozoic (Deng & Ding, 2015;Ding & Zhong, 2013;Royden et al., 2008), which strongly impacted the stability of environments in this region. The complex topography and highly variable climate in the EHHM are the main factors behind highly heterogeneous habitats or various niches observed in the region, which have been considered to account for the extremely rich biodiversity here (Wu, 1988;Xing & Ree, 2017;Zhang et al., 2009).
Generally, high mountains and wide rivers have significant effects on divergence and phylogeographic structure of species by creating strong geographic barriers for seed dispersal (Avise, 2000). Indeed, these effects have been demonstrated on some species in the EHHM region (Liu et al., 2009;Meng et al., 2007;Wen et al., 2014;Xu et al., 2010). However, some recent studies also have indicated that ecological factors, especially climate, might also play an important role in driving cryptic speciation or intra-specific divergence (Fan et al., 2013;Liu et al., 2013;Pinto-Carrasco et al., 2021;Yang et al., 2012;Zhang et al., 2016). Nevertheless, as Favre et al. (2015) reviewed, despite a growing number of studies, the specific mechanisms behind origin and evolution of biodiversity hotspots associated with the Tibetan Plateau remain poorly understood and further studies are needed.
For the patterns of biodiversity and climate change on the Qinghai Tibetan Plateau (QTP), the most serious dilemma we face is that we still cannot confirm when and how the QTP (the Himalaya Mountains (HM), one of the main constituent terranes of QTP), reached its modern elevation. There has been considerable diversity of opinions about the process and date of the plateau uplift (Deng & Ding, 2015;Renner, 2016;Royden et al., 2008;Spicer et al., 2021;Wang et al., 2014). These opinions often yield conflict inferences owing to different lines of evidences. Some researchers hold that the central region of the plateau rose to its present height as early as 40 million years ago (Ma), with subsequent outward extensions by the early Miocene (Rowley & Currie, 2006), while some other scholars suggested that the high central plateau was not formed until the Neogene (Su et al., 2019). Meanwhile, some other studies indicated that the uplift of the plateau reached over 4000 m average elevation only by 15 Ma (Coleman & Hodges, 1995;Spicer et al., 2003) or even 8-10 Ma (Deng & Ding, 2015;Harrison et al., 1992). In recent decades, drastically different hypothesis for recent and rapid uplift of the QTP at about 3.6 Ma was also proposed (Cui et al., 1996;Li et al., , 2015Shi, 1998). Thus, heated debates about the evolution history of QTP uplift still continue, in spite of unceasing accumulation of evidence from tectonics, fossils, isotopes, and climate simulations. Hence, in the region, studies of phylogeographic and phylogenetic patterns, which can link organic diversification with geological history and environmental change, are still needed.
Hippophaë is a genus of Elaeagnaceae and it may be one of the most appropriate taxa to clarify the mentioned above relationships.
The EHHM region is the diversification center of this genus, and six of seven species of the genus are distributed there (Jia & Bartish, 2018;Lian et al., 2006;. Hippophaë has a long evolutionary history (from the Eocene or early Oligocene) as indicated by both paleobotanic (Akkiraz et al., 2011;Miao et al., 2013) and molecular data (Bartish, 2016;Jia & Bartish, 2018). The range and the age of the genus both imply that the geological and climatic changes of the EHHM region should have influenced its evolution.
Our previous study on H. tibetana indicates a link between the uplift of QTP during the last 3~4 Ma, and that the climate change might play an important role in driving intraspecific diversification in H.
tibetana (Wang et al., 2010). Another species of this genus, H. gyantsensis (Rousi) Lian, may provide even better opportunity to test the presented above hypothesis. This species mainly occurs in the valleys of the middle Yarlung Zangbo River (YZR), only grows on riverbanks and floodplains, and ranges from 3000 to 4600 m in altitude (Lian et al., 2006). The distribution range of H. gyantsensis is now limited to the central Himalaya, not EHHM region, but we found several putative hybrid populations of this species in the Hengduan Mountains (personal communication), so we still focused on the EHHM region in this study. Age of the species was estimated be to the Early Miocene (Jia et al., 2016;Jia & Bartish, 2018) and its origin may be explained by several ancient hybridization events within the genus (Jia et al., 2016). Comparing with H. tibetana, this species has a continuous distribution and covers various climate zones from west to east (The Comprehensive Scientific Expedition to the Qinghai-Xizang Plateau, 1984). Importantly, there are no significant geographic barriers in its main distribution because its range is along the YZR valley, which is the channel for seed dispersal by birds. The fruits of Hippophaë plants are important food sources for many birds in the QTP (Lu et al., 2005) and thus can be dispersed over a long distance, which enables us to exclude the effects of strong Up to now, most phylogeographic studies in EHHM region concentrated on alpine plants, and very few of them studied the plants from the YZR there (Cheng et al., 2017;Wang et al., 2019). However, species growing in YZR are more sensitive to detecting the influence of climate dynamics. Here, through the inspection into the phylogeography of H. gyantsensis using both chloroplast DNA (cpDNA) and microsatellite fragments, we aimed to answer the following specific questions: (1) Does H. gyantsensis display significant phylogeographic structure? (2) Which role did climatic dynamics and the uplift of QTP play in shaping the phylogeographic structure of H.  Table 1.

| Chloroplast data analysis
The cpDNA trnT-trnF sequences were aligned with CLUSTAL X (Thompson et al., 1994), corrected manually in MEGA version 5 (Tamura et al., 2011), and then assigned to different haplotypes using DnaSP version 6 (Rozas et al., 2017) Jia and Bartish (2018). Simultaneously, we downloaded all the public trnL-trnF and ITS sequences (31 and 22, respectively) of H. gyantsensis from GenBank, of which 8 and 14 were respectively used in Jia et al. (2016). We used all these sequences for reconstructing the phylogenetic relationships within H. gyantsensis based on both trnL-trnF and ITS fragments. Phylogeny reconstruction with Maximum likelihood (ML) method was carried out using RAxML-VI-HPC (Stamatakis, 2006). Clade support was evaluated by bootstrap (1000 replicates). Bayesian Inference was conducted with BEAST version 1.7.4 (Drummond et al., 2012), employing the same model as used in ML analysis. The condition of Markov chain Monte Carlo (MCMC) was as follows: the total generation length was 10,000,000 generations, and trees were sampled each 1000 generations. The first 1,000,000 generations were discarded as burn-in. Convergences of each parameter were confirmed by the TRACER ver. 1.5, and ESS (effective sample size) of all parameters were larger than 200. The Network of cpDNA haplotypes was reconstructed with Network 10.3 (available at https://www.fluxu s-engin eering.com/share net.htm) using TA B L E 1 Geographic locations and genetic diversity of 21 populations of Hippophae gyantsensis based on chloroplast trnT-trnF region and 11 microsatellite loci, respectively. the Median-joining method (Bandelt et al., 1999) and MP calculation (Polzin & Daneshmand, 2003).
Haplotype diversity (H d ) and nucleotide diversity (π) of all the populations were estimated by DnaSP version 6. Average diversity within populations (H S ), total gene diversity (H T ), and two coefficients of population differentiation (G ST and N ST ) were calculated using PERMUT 1.0 (Pons & Petit, 1996). G ST (coefficient of genetic variation over all populations) is only based on allele frequencies while N ST takes the similarities of alleles into account, therefore significantly larger N ST than G ST suggests that similar alleles tend to be geographically closer, indicating significant phylogeographic structure. The significance was tested by 1000 permutations.
When the null hypothesis was not rejected, the time since expansion (t) was estimated according to the formula t = τ/2u (Rogers & Harpending, 1992), the value u = μkg, where μ is the substitution rate, k is the average sequence length, and g is the generation time.
Here, it was approximately estimated to be 5 years for H. gyantsensis (Bartish et al., 2006). In addition, neutral tests with Tajima's D (Tajima, 1989) and Fu's F S (Fu, 1997) were also conducted.
To define groups of populations, space analysis of molecular variance (SAMOVA) for cpDNA was conducted with SAMOVA 1.0 (Dupanloup et al., 2002). In this analysis, we set the number of groups (K value) from 2 to 6. According to the definition of groups by SAMOVA, an analysis of molecular variance (AMOVA) was also conducted using Arlequin 3.5 to calculate the molecular variance.
The significance tests were based on 1000 permutations.
We also tested whether the molecular clock was hold or not by using the BASEML program of PAML version 4 (Yang, 2007) on the basis of the ML tree topology. When the molecular clock hypothesis was not rejected, divergence times were estimated by the molecular clock. The substitution rate for the trnL-trnF region in Phylica (Rhamnaceae), 4.87 × 10 −10 substitutions per site per year (s/s/y), was adopted, which was summarized by Richardson et al. (2001) and used for H. tibetana (Wang et al., 2010). We used this calibration in our dating analyses. Jia and Bartish (2018) dated chloroplast phylogenies of Elaeagnaceae and Hippophaë by first resolving phylogenetic relationships within these taxa using a concatenated database of five chloroplast loci. They then calibrated these phylogenies by a fossil record of Shepherdia in North America from the Late Eocene.
According to these authors, mean age of the stem node of H. gyantsensis is about 20 Ma. We refrained from using this calibration because our interest was mainly in intraspecific differentiations within H. gyantsensis, which were found by Jia and Bartish (2018) to be much younger evolutionary events.

| Microsatellite data analysis
To measure the level of genetic diversity and genetic differentiation for microsatellite data, several programs were performed as and fixation Index (F) were estimated using GenAlEx v6.5 (Peakall & Smouse, 2012). In addition, we assessed population genetic structure under admixture model using the Bayesian method implemented in Structure 2.3.4 (Hubisz et al., 2009). The number of clusters (K) ranged from one to ten using 20 independent runs for each value of K. Each run comprised a burn-in of 10 6 generations, and followed by 10 6 MCMC steps. The optimal K was determined by log-likelihood value and ΔK statistics (Evanno et al., 2005). The population clusters were visualized using the software Distruct 1.1 ( Figure 2; Rosenberg, 2004). Subsequently, AMOVA was performed to partition total genetic variation within and among populations based on the result of structure analysis with 1000 permutations using Arlequin 3.5.
In addition, we evaluated the statistical support for eight alternative phylogeographic scenarios of the divergence history of H. gyantsensis ( Figure 3) based on the above-mentioned structure analysis using an Approximate Bayesian Computational (ABC) approach. We did not use the chloroplast data for ABC modeling due to insufficient variation in the sampled DNA fragment. Besides, unlike microsatellites, cpDNA represents only a small and specific part of the total genome of the species. Three groups were seta as western cluster (N 1 ), eastern cluster (N 2 ), and central cluster (N 3 ) after removing three mixed populations (i.e. P5, P7, and P11). The eight scenarios were designed in detail as below: in the first scenario, N 1 , N 2, and N 3 derived from an ancestral population at the same time (t 2 ). The other scenarios showed some (presumably most likely) of the possible patterns of divergence events among the three groups, of which the fourth scenario showed that N 3 was formed from a genetic admixture event between N 1 and N 2 at time t 1 . We gave a uniform prior probability and ran 8 × 10 6 simulations under each scenario using DIYABC 2.1.0 (Cornuet et al., 2014), of which 10% was used to estimate the relative posterior probability with 95% credible intervals via logistic regression and posterior parameter distribution (Appendix Figure S1;

| Climatic data analysis
To identify the climatic factors potentially associated with the divergence between groups of H. gyantsensis, we compared recent (c. 1950-2000) data of 19 BIOCLIM variables (Hijmans et al., 2005) between the groups. We extracted temperature and precipitation values from the BIOCLIM data sets with a grid size of 30″ (c.

| Microsatellite data
For the 11 microsatellite loci used here, genetic diversity indices were summarized for the population (   (Table 3). The differences of AMOVA result between microsatellite and chloroplast data were mainly resulted from their different evolutionary rate and history.
The ABC results suggested that a scenario that N 3 derived from historical admixture between N 1 and N 2 (Scenario 4) was the most probable model with 30.07% posterior probability (Figure 3).
According to the simulated results, there was no significant changes in effective population size between ancestral population and current populations. Posterior parameter estimates for Scenario 4 indicated that N 1 and N 2 diverged from each other at 65.5 ka (95% HPD: 17.5-132.0 ka), followed by a genetic admixture about 26.9 ka (Table 5).

| Climate heterogeneity
A total of 19 BIOCLIM variables were estimated in 21 populations and results were shown in Table 6. Out of the eight precipitation variables, annual precipitation (bio12), and precipitation seasonality (bio15) were significantly different between the two groups of H.
gyantsensis had a strong phylogeographic structure, and all sampled populations were divided into two distinct lineages, which occupy the western and eastern part of its range, respectively ( Figure 1b). In the revealed phylogeny, all the chlorotypes from the western group represented a strongly supported clade while all the chlorotypes from the eastern group formed another strongly supported clade, and the mean age of divergence between the two lineages was estimated at ~3.6 Ma (Figure 1a). This estimate is within the range (1.2-3.9 Ma) reported by Jia and Bartish (2018) for the crown node of the species. Therefore, both age estimates (the earlier study and ours) placed the earliest diversification in  Appendix Figure S4) These results, a long temporal divergence between the lineages, a short spatial distance, and obvious opportunities for long-distance dispersal of seeds between the ranges of two lineages, indicate that geographic barriers are unlikely to drive the formation of phylogeographic structure of this species and the divergence of the two lineages.

| Effects of climate factors on the formation of phylogeographic structure of H. gyantsensis
Previous studies have found significant phylogeographic structure in plant species from the EHHM region (Ge et al., 2005;Li et al., 2011;Qiu et al., 2011). The formation of such structure is often attributed to vicariance, as gene flow is often blocked by geographical barriers (e.g. high mountains). However, some many studies also found that the split between lineages does not strictly follow the geographic barriers, but geological isolation and ecological factors, especially climate, together promote the divergence of several species from this region (Fan et al., 2013;Liu et al., 2013;Yang et al., 2012;Zhang et al., 2020). In the study on H. tibetana, the split between lineages coincided with 400 and 600 mm annual precipitation lines, also suggesting that climate played an important role in driving intraspecific diversification in H. tibetana (Wang et al., 2010 Abbreviations: N a , ancestral population size; N 1 , N 2 , and N 3 represent effective population size of western cluster, eastern cluster, and central cluster, respectively. t, divergence time or admixture time of different clusters; μ, the mutation rate; r a , admixture rate. TA B L E 6 BIOCLIM variables in H. gyantsensis populations and the results of two-tailed t-tests between groups and analysis of variance (ANOVA).  (Widrlechner, 1997). Figure 5 shows plant hardiness zones in the EHHM according to China Plant Hardiness Zones (CPHZ) published by Widrlechner (1997). It was found that almost all populations of the eastern group were located in CPHZ 7 (−17.7 to −12.3°C) while those of the western group in CPHZ 5 or 6.
In the QTP, precipitation has also been considered as a key climate factor to determine the divergence of species, and annual precipitation lines of 400 and 600 mm were two important boundaries for different lineages of H. tibetana (Wang et al., 2010).
Between the western and eastern parts of H. gyantsensis, precipitation seasonality (bio15) have also showed significant differences in average values (Table 6), and their boundary was very near the annual precipitation line of 400 mm (Figure 4). In this study, the relevance of genetic divergence of H. gyantsensis with June precipitation was also analyzed due to its potential importance and for obtaining more accurate data. As the beginning of the India monsoon every year, June precipitation is not only very closely related to the precipitation seasonality on the QTP, but is also an important factor for the growth and development of H. gyantsensis since this plant is flowering in this month (Crimmins et al., 2011).  Widrlechner (1997).

| Recent regional uplift of the QTP and effects on the formation of phylogeographic structure of H. gyantsensis
As discussed before, the climate heterogeneity in the reaches of YZR is the most important and direct factor to drive the divergence and phylogeographic structure of H. gyantsensis. This heterogeneity is mainly determined by two key factors: the HM and Indian monsoon (IM). The climatic heterogeneity in the YZR occurred only when the HM reached a particular elevation and IM was driven up to certain strength. Therefore, the origin of the above heterogeneous climate must be accounted for timing of the rise of the HM, and the establishment of the IM.
Although both the time and course of their uplifts have been controversial (Deng & Ding, 2015), most of available so far evidences supported the suggestion that QTP and the HM have different histories (Deng & Ding, 2015;Li et al., 2015;Spicer et al., 2021;Su et al., 2019;Wang et al., 2014). The QTP might had risen to the recent height as early as 40 Ma (Tada et al., 2016;Wang et al., 2014;Wang, Schluetz, & Liu, 2008;Wang, Zhao, et al., 2008), but the uplift of the HM was s likely later than that of the QTP (Deng & Ding, 2015;Ding et al., 2017). Some studies on the basis of oxygen isotope indicate that the HM had reached its modern height or even higher than present by 20 Ma (Ding et al., 2017;Gébelin et al., 2013). However, some other studies suggested that before ~10 Ma, the HM were still dissected by some of the deepest and most impressive gorges on Earth and rivers rising in Gangdese (Transhimalaya) would have flowed across the HM and flowed south into Indian plain (Tremblay et al., 2015). In addition, fossils in Gyirong and Zhada showed these regions (the north and center of the HM) were still warm and moist in the late Miocene and even in the early Pliocene (Huang et al., 2020;Huntington et al., 2015;Ji et al., 1980), strongly suggesting the HM was unlikely to reach the recent height before the late Miocene.
Furthermore, the massive conglomerate which has long been regarded as the product of orogeny around the QTP margin began to deposit at about 3.6 Ma, implying the central HM underwent a rapid rise in the Pliocene (Li et al., 2015). All these studies showed that different evidence from different regions of the QTP may reflect discrepant evolution history. This interpretation makes sense if: (1) the QTP is not a single geological entity but a fusion of several accreted terranes, and they evolve in piece-meal manner (Fielding, 1996;Spicer et al., 2021;Su et al., 2019;Tapponnier et al., 2001;Wang, Schluetz, & Liu, 2008;Wang, Zhao, et al., 2008); (2) different types of data measure different aspects of the topography. For example, fossils are inclined to reflect low elevations, while isotopes tend to reflect high altitude (Botsyun et al., 2016;Spicer et al., 2021). Studies in phylogeography (microevolution) and phylogeny (macroevolution) frequently relate divergence-time estimates to paleogeographic or climate events, often infer the abiotic factors making for clade diversification or speciation (Favre et al., 2015;Renner, 2016). Most phylogeographic, phylogenetic and biogeographic studies on taxa of the QTP were linked to rapid and recent uplift of the QTP (the late Miocene and later; e.g., Cheng et al., 2017;Meng et al., 2017;Wang et al., 2005Wang et al., , 2010 (Deng et al., 2011;Wang et al., 2012;Zhang et al., 2012;Zheng et al., 2000). These results indicated that the possibility of species differentiation, especially those with current ranges at high altitudes and adapted to alpine ecosystems, was affected by the younger geological events.
In addition, the IM is another key factor to determine the climate of the EHHM, which has been considered to be closely related to the uplift of the HM and QTP (Boos & Kuang, 2010;Tada et al., 2016).
Regarding the onset of IM, more recent studies suggested that this appears to have begun during the late Middle Miocene (~12.9 Ma; Betzler et al., 2016) or even earlier (Guo et al., 2002;Huber & Goldner, 2012) and summer monsoon was in its full strength in the late Miocene (~7 Ma). Despite some controversy in understanding of the relationship between IM development and HM uplift, the perception that IM intensification occurred in the Late Pliocene is rarely disputed Lu et al., 2020;Zhang & Liu, 2010;Zheng et al., 2004). By this time, the divergence of the two lineages of H. gyantsensis may had occurred, according to our estimates.
This case adds therefore support to the hypothesis of partially recent and rapid QTP uplift. The concordances between geological, climatic, and biotic processes imply that the likely recent uplift of HM acted as creator of a high climatic heterogeneity in this region.
The reinforcement of Indian monsoon can also be associated with the divergence of this species. Moreover, some species divergences could be better explained by even more recent geological uplift (Wang et al., 2010;Xing & Ree, 2017 (Li et al., 2015), further indicated strong links between geo-climatic and biotic processes in the region.

| Demographic history and population genetic admixture
In this study, we simultaneously used chloroplast and microsatellite data to estimate the demographic history of H. gyantsensis.
Theoretically, these two types of data can respectively explain ancient and more recent evolutionary history due to the differences in mutation rate (Bai & Zhang, 2014). Besides, nuclear and chloroplast genomes are transmitted between generations in different ways.
Only nuclear genome information is carried by pollen, while both nuclear and chloroplast genomes are dispersed via seeds in Hippophaë , as in most other angiosperms. The difference in transmission routes from parent to offspring generations for different genomes (chloroplast and nuclear) can lead to differences in effective population sizes between exclusively maternally (via seeds) transmitted chloroplast genome, and both maternally and paternally (via pollen) transmitted nuclear genome. For strictly outcrossing dioecious Hippophaë, this means, in theory, that (all else being equal) effective population sizes measured from sequences of chloroplast genes are one-quarter of correspondent values based on nuclear genes (Birky et al., 1983).
For chloroplast data, we found east-west phylogeographic break in this species, and identified population expansion in the eastern cluster at about 0.12 Ma, when the TP was in the last interglacial period (MIS 5, Cui et al., 2011;Zheng et al., 2002). Similarly, a strong population genetic structure was also estimated from microsatellite data, but with three clusters (i.e. eastern cluster, central cluster, and western cluster). However, these three clusters have not displayed signs of demographic expansions in recent time (Table 5).
Furthermore, ABC modeling indicated that the origin of central cluster (P8-P10, P11-P12) was a result of admixture between clusters eastern and western at 27 ka (Table 5) (Cui et al., 2011). This westward dispersal event can be evidenced by the genetic composition of P7 and P11. These populations possessed chlorotypes of the western group, while microsatellite data revealed genetic components of both eastern and western clusters. The significant phylogeographic break in cpDNA data but not in nuclear data for most species was often explained by genetic introgression across spatially narrow admixture cline (Barton & Gale, 1993;Cheng et al., 2021;Li et al., 2019). The indicated that demographic expansion of eastern cluster of H. gyantsensis likely provided opportunity for the genetic introgression of mixed populations. Additionally, the time estimates of population expansion and genetic admixture event in this study were both traced back to warm inter-glaciation during the Quaternary glaciation in China. Taken together, these results are potentially highlighting the importance of consequences of Quaternary climatic fluctuations for regional biota (i.e. glacial and interglacial isolation) to the recent evolutionary history of H. gyantsensis (Liu et al., 2014;Wen et al., 2014).

| CON CLUS ION
Our chloroplast and microsatellite analyses of H. gyantsensis clearly revealed that range of H. gyantsensis was most likely limited to the middle YZR in the Late Neogene. We also found significant phylogeographic structure in the species with strong spatial differentiation into eastern and western groups. Age of divergence between the two main lineages of the species could be traced to ~3.6 million years ago. A combination of the recent regional uplift of the China (Grant Nos. 31860167, 91131901, 41061007, 31160046, and 31270407).

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.